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ABSTRACT 

Using two- and three-dimensional hydromagnetic simulations for a range of different flows, 
including laminar and turbulent ones, it is shown that solutions expressing the field in terms of 
Euler potentials (EP) are in general incorrect if the EP are evolved with an artificial diffusion 
term. In three dimensions, standard methods using the magnetic vector potential are found 
to permit dynamo action when the EP give decaying solutions. With an imposed field, the 
EP method yields excessive power at small scales. This effect is more exaggerated in the 
dynamic case, suggesting an unrealistically reduced feedback from the Lorentz force. The EP 
approach agrees with standard methods only at early times when magnetic diffusivity did not 
have time to act. It is demonstrated that the usage of EP with even a small artificial magnetic 
diffusivity does not converge to a proper solution of hydromagnetic turbulence. The source of 
this disagreement is not connected with magnetic helicity or the three-dimensionality of the 
magnetic field, but is simply due to the fact that the nonlinear representation of the magnetic 
field in terms of EP that depend on the same coordinates is incompatible with the linear 
diffusion operator in the induction equation. 
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1 INTRODUCTION 

In the past few decades magnetic fields have become an inte- 
gral part of many branches of observational and theoretical astro- 
physics. This is because in virtually all astrophysical bodies the 
electrical conductivity is large enough to support electric currents 
and hence magnetic fields. Furthermore, virtually all astrophysi- 
cal flows produce dynamo action allowing part of the kinetic en- 
ergy to be channelled throu gh the magnetic energy reservoir be fore 
it is being dissipated (see iBrandenburg & SubramaniarJ 120051 for 
a review). Simulating such flows on the computer can become a 
serious challenge, in particular if one wants to reach large mag- 
netic Reynolds numbers and if one wants to represent huge den- 
sity contrasts that are typical for self-gravitating centrifugally sup- 
ported structures such as galaxies. The same challenge is met in 
cosmological simulations that describe the formation of galaxy 
clusters and even the formation of galaxies. Many such simula- 
tions have been pe rformed using smoothed particle hydrodynamics 
dDolag et alj|2002k Its Lagrangian nature is well suited for han- 
dling self-gravity dMonagha n 1992). However, incorporating mag- 
netic fields into such simulations has proved challenging. A pos- 
sible solution to this problem may be the use of Euler potentials 
jPrice & Bat3l2007l : iRosswog & PricefcOOl 120081) . 

On a number of occasions the use of Euler potentials 
(EP) ha s proved us e ful in astrophysics and magnetohydrody- 
namics dSweetl 1 19501: iDungeyl 1 19581 ; ISternl 1 197(A ISakurail [19791 ; 
lYahalom & Lvnden-B ell 2006). In this approach the magnetic field 
is written as 



Vq x V/3, 



(1) 



where a and (3 are the EP. Until recently, the use of EP has only 
been modestly popular, because the nonlinearity of such a represen- 
tation of B can lead to difficulties in repres enting a r bitrary initial 
conditions. Furthermore, as pointed out by Moffat3 dl978f) . mag- 
netic fields with linked or knotted B lines cannot be represented 
with single-valued differentiable EP. Another problem is that one 
has the evolution equations for a and [3 only in the ideal case, i.e. 
when the resistivity vanishes. In that case one has to solve just two 
simple advection equations, 



Da/Dt = and D/3/Dt = 0. 



(2) 



Here,D/Df = d/dt+U-V is the advective derivative and U is the 
velocity. In recent years the use of Euler potentials has become in- 
creasingly popular in SPH simulations, because the evolution equa- 
tions for a and /3 imply that the values of a and /3 are simply kept 
fixed at all times. Several tests have suggested that the use of EP 
can be superior to solving for B because of t he difficulty in pre- 
servin g V • B — in numerical simulations dDolag & Stasvszvnl 
2008). Differences between the two results could therefore readily 
be explained in terms of V B not being zero in the latter approach. 
However, this does not eliminate concerns about the correctness of 
solutions obtained with the EP method compared to other methods 
that also preserve V ■ B = 0. One such method is to solve for the 
magnetic vector potential (A method). In this paper we compare the 
two methods for a range of different flows. 
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2 EULER POTENTIALS IN SIMULATIONS 

Kota rba et all feOOSl, discussed the fact that the magnetic helicity 
vanishes in the EP representation, and so it is clear that this method 
is not well suited for studying helical or a effect dynamos that tend 
to produce magnetic fields with finite magnetic helicity. However, 
we still do not know what the magnetic field will be in such a case, 
and whether the EP method can still be useful for studying other 
types of dynamos, or at least other types of turbulent magnetohy- 
drodynamic (MHD) flows. 

The goal of this paper is to compare the evolution of the mag- 
netic field in simulations using the EP method on the one hand 
and the magnetic vector potential method (A method) on the other. 
We do this by solving the equations for both methods at the same 
time. Of course, in the majority of cases, sharp gradients may de- 
velop eventually. This is when numerical methods for solving the 
equations of ideal MHD break down. It has then been customary 
to include artificial diffusion in the evolution eq uations for a and 
(3, i.e . one c onsiders solutions of the equations (Rosswog & Price 
l2007ll2008l) 

D/3 



Da ^2 



Dt 



(3) 



where r\ is the magnetic diffusivity. In the two-dimensional case 
with B = B(x, y) and B z = we can write B — V x (A z z), 
where z is the unit vector in the z direction, and A z obeys the 
uncurled induction equation, which can be written as 

DA 



Dt 



71^ A z 



(4) 



To compare with the EP method, we choose /3 = z and write 
A = aV/3 = az, where z is the unit vector in the z direction, 
so we have A z — a(x, y, t), and thus the evolution equation for a 
becomes identical to that for A z , even when r\ 7^ 0. One can also 
write A = — /3Va, which agrees with the previous formulation 
after adding the gradient of a/3, which does not affect the B field. 

In order to facilitate direct comparison between the EP and 
A methods, we solve numerically Equation (O together with the 
equation for the A method (Appendix lAl). 

M = - A -(VUf+r,V 2 A, (5) 

where we have assumed 77 = const. We emphasize that the velocity 
U enters Equations (3]l and Q also through the D/Dt derivative, 
and that the equations for both approaches are equivalent in the 
special case of 77 = 0. Indeed, if we insert a symmetrized represen- 
tation, 



A = i(aV/3-/3Va) 
into Equation {5}, we obtain 



(6) 



("57 V/3_ (,1m Vn ~ J? + v " <7 > 



/D/3 



where R stands for a residual term, and (f> is given by 

</> = §(a/3 - /3d) + i//(aV 2 /3 - (3V 2 a) -UA. 



(8) 



The R term vanishes for 77 = 0, but is finite with magnetic diffu- 
sivity and is then given by 



R = ij(Va ■ V)V/3 - r?(V/3 • V)Va. 



(9) 



Note that the term can be removed from Equation ([7} by a 
gauge transformation; see Appendix iBlfor the derivation. However, 
the R term cannot be removed and, moreover, it has the same high- 
est order of derivatives as the terms on the LHS of Equation Q, 



so R is in general not small. This is exactly the reason why the in- 
troduction of artificial diffusion is in general not permissible. The 
hope is, of course, that in the limit 77—5-0 the EP and A methods 
give still reasonably similar results. In order to illustrate when this 
is the case, we consider in the following different flow fields. 



3 CHOICE OF FLOW FIELDS 

We first consider the case where U is a given function and turn 
then to the case where U is obtained by solving the momentum 
and continuity equations. In the former case we restrict ourselves 
to flows of the form 



U = V x ipz + <f>z, 



(10) 



where ip — ?p(x,y,t) and (f> = cj>(x,y,t) are prescribed func- 
tions that will be defined below. In the latter case we consider the 
compressible equations with an isothermal equation of state, so the 
density p is proportional to the pressure, which is then given by 
p — pc 2 , where c s = const is the isothermal speed of sound. The 
governing equations are then 



Dlnp 

Dt 
DU 



= -v-v, 



Dt 



= -c 2 Vlnp + f + F vis 



(11) 



(12) 



where f vise = P~ X V ■ 2pvS is the viscous force, v is the kine- 
matic viscosity, Sij = \{U%,j + Uj,i) — ^Sij'V ■ U is the traceless 
rate of strain tensor, and / is a nonhelical random forcing function 
consisting of plane transversal waves with random wavevectors k 
suc h that [fcj lies in a band around a given forcing wavenumber 
k { dHaugen et al.ll2004l) . The vector k changes randomly from one 
timestep to the next. The forcing amplitude is chosen such that the 
Mach number Ma = it rms /c s is about 0.1. 

The total system of equations consists of Equations I0 and (O 
together with Equations i ll It and H21 . In all cases the magnetic 
field is considered infinitesimally weak, so that the Lorentz force 
can be neglected. These equations were solved using the PENCIL 
which is a high-order finite-difference code (sixth order in 
space and third order in time) for solving the compressible MHD 
equations. The code came with a routine that solves two passive 
advection-diffusion equations that were invoked by compiling with 
CHIRAL=chiral, which is a routine that was originally designed 
for another purpose to des cribe the spontaneous chiral sym metry 
breaking in biomolecules (Brandenburg & Multamaki 2004). Ad- 
ditional diagnostics for monitoring the magnetic field and the cur- 
rent density have been added to this module for the purpose of this 
paper. 

Initial conditions are generated by setting first a and j3, and 
then calculating A from Equation l[6}. We consider cubic domains 
of size L 3 using triply-periodic boundary conditions in all cases, 
except the first one which is a two-dimensional case where we as- 
sume perfect conductor boundary conditions. In either case, the 
magnetic helicity, H = J A ■ B dV, is gauge invariant, i.e. the 
transformation A — > A' + VA does not change the value of H. 
This is because 



VK-BAV = js AB AS- I XV Bd\' 



(13) 



1 http ://pencil-code . googlecode.com/ 
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vanishes owing to the condition V ■ B = 0, and there is no surface 
term for periodic domains or perfectly conducting boundaries. So, 
the statem ent that in th e EP approach A ■ B — is merely a gauge 
condition dSternll970i) does not change the fact that we always have 
H — 0. On the other hand, the current helicity, J J ■ B dV, can 
well take values different from z ero, as has bee n utilized in the cal- 
culation of force-free equilibria (Sakurai 1979). In the A approach 
H can generally be different from zero. However, if H = ini- 
tially, then H can become different from zero only through resis- 
tivity. This is a consequence of periodic or perfectly conducting 
boundaries. 




4 RESULTS 

4.1 Wind-up by a two-dimensional eddy 

We consider first the wind-up of an initially uniform magnetic field, 
B — Box, or a — Boy and f3 — z in the EP formulation. We 
choose a flow with a single eddy given by Equation ( IIOK i.e. 



ip(r) = (Uo/k) cos kr, </>(r) = eUoip{r), 



(14) 



where r 2 — x 2 + y 2 in a domain —L/2 < x, y < L/2 and 
k = tt/L. For e = 0, this flow was use d earlier to compute the 
magn etic field evolution in ideal MHD ferandenbur g & Zweibell 
1994), so we were able to compare our results with theirs in the 
ideal case. 

We adopt perfect conductor boundary conditions, which corre- 
sponds to keeping the values of a and /3 on the boundaries equal to 
their initial values. However, it is advantageous to subtract out the 
linear gradients of a — cto + oti and /3 = (3 + Pi and solve only 
for the departures ai and /3i, whose values vanish on the bound- 
aries. In our case the imposed gradient fields are ao = Boy and 
f3o = 2, so the relevant evolution equations are 



Dai TJ „ , ^2 W 



- = -U z +rf7 2 Pi. 



(15) 



The result is shown in Figure Q] for the ideal case, r\ = 0, with 
e = 2/tt and different resolution. One sees clearly that B Ilas in- 
creases linearly with time while the current density J — V X B/y,o 
(with fio being the vacuum permeability) increases quadratically 
with time. In B rms the differences between EP and A methods are 
small, which is why we plot in the second panel the maximum value 
of | J\. Departures from the more accurate solutions obtained at the 
next higher resolution appears roughly at the same times, but are 
seen more clearly in J m ax than in B Tms . The linear and quadratic 
scalings for B and J, respectively, are well reproduced by either 
method provided the resolution suffices to resolve the progressively 
finer structures as time goes on. Looking at the plot of J max , one 
can conclude that one may need slightly more points with the EP 
method than with the A method. 

For e = the EP method gives correct results even in the case 
of finite magnetic diffusion, as expected based on the equivalence 
of the underlying equations in that case. This is connected with the 
fact that the flow is two-dimensional and confined to the plane only, 
However, when Owe have U z (x,y) ^ and B z (x,y) 7^ 0, 
and hence f3i 7^ 0. In this case, the R term is in general non- 
vanishing, and so Equations <[3j and l |15t are then no longer equiva- 
lent to Equation l[5}, even though the flow and the field depend only 
on two spatial coordinates. This is demonstrated in Figure|2] where 
we plot the time dependence of the current helicity, {J ■ B), in runs 
with zero and finite values of rj. Note the mutual departure of the 
two methods after some time when 77 7^ 0. 
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Figure 1. Evolution of B r ms/Bo (upper panel) and (Jmax/Jo) 1//2 (lower 
panel) for different resolutions with r\ = and e = 2/tt. Here, Jo = 
kBo/fio has been used for normalization. In both plots dashed lines gives 
the ideal scalings, i.e. linear for B and quadratic for J. 
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Figure 2. Evolution of the current helicity (J ■ B) for e = 2/tt with 
r) = (where A and EP methods both give the same quadratic scaling; 
dashed line) and r\ 7^ (where the two methods disagree after some time). 
Again, Jo = kBo / fio has been used for normalization. 



These experiments have demonstrated that our implementa- 
tion of the EP method along with the corresponding diagnostic 
tools give agreement with the A method, even when e / 0, pro- 
vided r\ — 0. However, for r\ 7^ the two methods only agree 
when e — 0. It may appear that the disagreement is connected with 
the occurrence of current helicity, but this is not the case, as will be 
discussed at the end of the paper. 
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Figure 3. Comparison of the evolution of B rms in a Roberts flow for meth- 
ods A and EP for a smooth initial condition (upper panel) and a random one 
(lower panel) for resolution of 128^ meshpoints and i? m = Uo/r)k = 200. 
Both plots are double-logarithmic, so as to see more clearly the mutual de- 
partures of the two solutions at early times. The insets give the more usual 
linear-logarithmic representation showing clearly the exponential growth 
of the A solution at later times. The dash-dotted line departing from the EP 
line at t = 3/Uok is the result of the A method, but with an initial condi- 
tion calculated from the EP solution at that time (Run A2), as opposed to 
the initial time (Run Al). 

4.2 Roberts flow dynamo 

Next we discuss the iRobertj d 19721) flow given by Equation d 1 01 > 
with 

t/j(x,y) = (U /k) coskxcosky, <j>(x, y) = k f i[)(x, y), (16) 

in the domain — L/2 < x,y < L/2, with k — 2n/L and kf = 
\p2k. The Roberts flow is one of the simplest flows that produce 
dynamo action. The dynamo is however a slow one, i.e. its growth 
rate goes to zero in the limit rj — » 0. The critical value of rj is 
?7crit = 0.181f/o/fc, so the critical value of the magnetic Reynolds 
number is i? m = Uok/ij cr i t = 5.52. 

We have considered two different initial conditions, a smooth 
one given by a = cos ky and /3 = cos kz, and a random one 
where a and j3 are given by independent random functions. The 
results are shown in Figure [3] For smooth initial fields the EP and 
A methods agree up to 8 time units (Wok = 8). This suggests that 
the EP method gives valid results only when magnetic diffusion did 
not yet have time to act. The A method shows that dynamo action 
commences after 100 time units, while the EP method gives only 
decaying solutions. 

For random initial fields dynamo action occurs earlier, after 
about 10 time units, but the growth rate is the same as for smooth 
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Figure 4. Comparison of the evolution of B rlnB for the modified Galloway- 
Proctor flow with point-wise zero helicity for methods A and EP using 256 3 
meshpoints and i? m = 10 4 . Note the power law scaling for the EP method 
and the exponential scaling for the A method. 

initial conditions. The reason for the difference in the onset of the 
exponential growth is that the eigenfunction of the dynamo mode 
overlaps poorly with the smooth initial condition, and does better so 
with random initial conditions. However, for random initial fields 
(with a spatially white noise power spectrum) the EP and A meth- 
ods give different results from the very beginning (Runs Al and 
EP). This is because of large discretization errors associated with 
the numerically different representations of white noise spectra. In 
order to check this we have calculated a new initial condition from 
the EP solution at time t = 3/Uok, when the field has become 
sufficiently smooth to be accurately represented by both methods. 
Now there is initial agreement, but it is still followed by a departure 
immediately afterwards (Run A2). 

The results demonstrate quite clearly the difference with the 
EP meth od in handling helica l dynamos, just as anticipated previ- 
ously by iKotarba etaij 12009). However, it is still unclear whether 
the helical Roberts dynamo is just an exception, or whether the dif- 
ferences are of more general nature. 

4.3 Flows with point-wise zero helicity 

The problem with the Roberts flow is two-fold. Firstly, it is clear 
that the dynamo produces a large-scale field of Beltrami type and 
is therefore helical. This is impossible to represent in terms of EP. 
Secondly, the dynamo does not exist in the limit ij — » 0, which is 
the only case where there is hope that the EP method can work. 
The latter problem could potentially be alleviated by choosing a 
flow that permits fast dynamos, where the growth rate remains finite 
in the limit rj — ► 0. However, this may not be true if 77 — > is 
a singular limit, w hich is different from the case r\ = 0. Time- 
dependent flows of iGallowav & Proctor] i 19921) type tend to be fast 
dynamos. An example of such a flow that has also point-wise zero 
kinetic helicity is given by (see lHughes et al . 1996) 

i>{x,y,t) = y/s/2(U Q /k)[coskX(x,t) +sinkY(y,t)], (17) 

4>{x, y,t) = k sin kX(x, t) cos kY(y, t), (18) 

kX — kx + cosut, kY — ky + sinwi, (19) 

in the domain —L/2 < x, y < L/2, with k — 2n/L. 

In Figure [4] we show an example for R m — Uo/rjk = 10 4 . 
Again, it turns out that the EP method does not give solutions that 
are compatible with those of the A method. It turns that, while for 
the A method the field grows exponentially like B IU1B ~ e xt , for 
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Figure 6. Comparison of the evolution of B rms in nonhelical turbulence 
for methods A and EP using 128 3 meshpoints at i? m = 80 (Runs Al and 
EP1) as well as i? m = 160 (Runs A2 and EP2). Note that the growth rate 
for Run A2 is slightly larger than that for Al, while the decay rates for EP1 
and EP2 are the same. 
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dependence of the exponents A and a characterizing the evo- 
At for the A method and S rms ~ t" for the EP method 



Figure 5. R 

lution of firms 

for the modified Galloway-Proctor flow with point-wise zero helicity for 
methods A and EP using 256 3 meshpoints. 



the EP method the field decays algebraically like B rms ~ t~ a . In 
Figure [5] we plot the dependence of A and a on i? m . It turns out 
that A seems to converge to a finite value (for the A method), and 
so does a (for the EP method), confirming that the functional forms 
of the time dependencies for the A and EP methods are indeed dif- 
ferent even for large values of i? m . 



4.4 Nonhelically forced isotropic turbulence 



The flows considered in Sections [4. 2l and l431 are laminar. Another 
example of fast dynamo action, where the growth rate is compara- 
ble to the inverse turnover time, is isotropic turbulence. This is also 
the example that is closest to the application to turbulence in galaxy 
clusters. In that case there might be a chance to see a tendency to- 
ward dynamo action with the EP method when r/ — > 0. Unlike 
dynamos with helicity, we can only expect the magnetic field to 
have length scales smaller than the energy-carrying scale. 

We consider here the case kf/ki = 3, with fci = 2n/L 
and a magnetic Reynolds number _R m = u ims /r]kf « 80, using 
v = T). The result is shown in Figure [6] lust like in all previous 
cases, there is a stark difference in the evolution of the magnetic 
field computed with the A and EP methods. With the A method 
we reproduce ex ponential growth cons i stent with earlier findings 
in the literature Jcho & Vishniadl200fi ISchekochihin et alj|2002t 
lHaugen et al.l2003l) . while the EP method gives results that bear no 
resemblance with those where small-scale dynamo action is possi- 
ble. The same is true of cross-sections of the field; see Figure [7] 
This strongly suggests that the EP method does not provide a solu- 
tion that is close to the expected one, except for the case of a planar 
flow that depends only on two coordinates. 



4.5 Spurious growth 

In the early days of dynamo theory there have been cases of grow- 
ing solutions that later turned out to be spurious due to lack of 
resolution. To demonstrate this in the present case, we present in 




EP method 
WD 




Figure 7. Comparison of cross-sections of B z (x, y) for methods A and 
EP from Runs Al and EP1 after 200 time units, using 128 3 meshpoints at 
Rm = 80. Light (yellow) shades indicate positive values and dark (blue) 
shades indicate negative values. Note the absence of any resemblance be- 
tween the two fields. 



Figure[8]a solution with f] = 0, keeping the fluid Reynolds number 
equal to 80, as in Figure [6] 

The A and EP solutions show obvious signs of insufficient 
resolution with oscillation on the scale of the mesh. Nevertheless, 
both solutions show exponential growth with the same growth rate, 
which is spurious given the presence of oscillation on the mesh 
scale. This illustrates the importance of considering the dependence 
of the solutions on r/, as was done in Section |4~4] In that case it 
turned out that the growth rate increases slightly with R m , but this 
behavior was not reproduced by the EP method. 



4.6 MHD turbulence with imposed field 



The problems considered in Sections |4. 2H4. 5 1 had to do with dy- 
namo action. This raises the question whether discrepancies be- 
tween the A and EP methods also exist in other cases where there 
is no dynamo action. As an example we now consider nonhelical 
turbulence in the presence of an imposed field using ao = Boy 
and Po = Z as initial fields, similar what was done in Section |4~T| 
The energy density of the imposed field Bo is comparable to the 
kinetic energy density. This is strong enough to dominate over dy- 
namo action and may even suppress it. Here we choose the forcing 
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Figure 8. Spurious growth of B rms to the right of the vertical line 
(titrmsfcf > 25.5), for nonhelical turbulence and methods A and EP using 
128 3 meshpoints with no resistivity (77 = 0), and a fluid Reynolds number 
of 80. 



wavenumber to be kt/ki = 1.5. The fluid and magnetic Reynolds 
numbers are again around 80. 

We consider both the kinematic case without feedback onto 
the flow and the dynamic case where the Lorentz force per unit 
mass, J x B / p, has been added to the rhs of Equation J 1 2b sepa- 
rately for the A and EP methods. The results are shown in Figure[9] 
where we plot magnetic power spectra for the two methods. It turns 
out that with the EP method, both the kinematic and dynamic cases 
yield excessive spectral magnetic energy at smaller scales (larger 
wavenumbers) compared to what the A method gives. Note also 
that with the EP method the resulting Lorentz force is weaker than 
with the A method, making the discrepancy even more pronounced 
in the dynamic case. With the A method kinetic and magnetic en- 
ergy spectra are in approximate equipartition with each other, while 
with the EP method the magnetic field exceeds the spectral kinetic 
energy at progressively smaller scale. 



5 DISCUSSION 

Having demonstrated that under a number of circumstances of 
practical interest the EP method is unable to provide a meaning- 
ful trend when artificial diffusion is added, one wonders whether 
the source of this failure can be identified more precisely. In par- 
ticular, we want to know whether this failure is connected with the 
inability to represent magnetic fields with helicity, or whether it is 
connected with the fact that the magnetic field is three-dimensional. 

In order to address these issues, we consider now a simple 
decay problem with U — and look for solutions of the equations 

OA ^ 2 . dB 
that disagree with solutions of the equations 



(20) 



da _ 2 d(3 _ 2 „ 



(21) 



even though the initial conditions obey B = Va x V/3. The 
essence of the problem can already be demonstrated with a non- 
helical field in two dimensions. An example is 



a = — cos ky, f3 = cos kx sin ky, 
which gives 

B(x,0) = (0, 0, k 2 sin kx sin 2 ky) 



(22) 



(23) 
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Figure 9. Magnetic power spectra for the kinematic case (upper panel) and 
the dynamic case (lower panel) for the A method (solid lines) and the EP 
method (dashed lines). In the dynamic case the kinetic energy spectra are 
also shown as thin solid and dashed lines for the A and EP methods. Here, 
kf/ki = 1.5 and R m = 80. 



as initial field. Note that 

poJ{x, 0) = (2 sin kx cos ky, — cos kx sin ky, 0)k sin ky, (24) 

so J ■ B = 0. With periodic boundary conditions, Equation J2U 
results in exponential decay of a and j3 while, owing to the non- 
linear representation of B — Va x V/3, the B field shows a 
non-exponential decay; see Figure [10] This problem is also clear 
from the fact that the R term in Equation does not vanish. This 
is generally a consequence of a and f3 being simultaneously depen- 
dent on the same coordinates (in this case both a and (3 depend on 
y). Alternatively, if we choose a — a(y) and f3 = j3(x) with 

a = ifci/ — i sin 2ky, (3 = cos kx, (25) 

which also results in B(x, 0) given by Equation l |23| l, then R = 
and a(y, t) shows a non-exponential decay — compatible with the 
correct solution of B(x, t). 

In general, a and f3 are functions of all three coordinates, so 
the R term in Equation l[9} does not vanish and the EP method with 
artificial diffusion will give wrong results. Thus, we can say that the 
failure of the EP method in the presence of artificial diffusion is not 
related to magnetic helicity nor to the three-dimensionality of the 
magnetic field, but simply to the fact that the nonlinear representa- 
tion of the magnetic field in terms of independent functions a and 
j3 is incompatible with the linear diffusion operator. 



6 CONCLUSIONS 

The EP method can give reliable results when r\ — 0, i.e. when 
one is interested in solutions to the ideal MHD equations. How- 
ever, in practice this is not possible, especially when the flows are 
turbulent, because energy must be dissipated at the smallest scale. 
When one allows for magnetic diffusion to be present, there is no 
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Figure 10. Decay of B ImB for solutions of Equations (20} and )2U using 
as initial conditions those given by Equations (22) and (23), respectively. 



agreement between the EP and A methods. As a consequence, it is 
impossible to use the EP method to study dynamos. Even fast dy- 
namos, which have finite growth rate in the limit 77 — 5- 0, cannot be 
modelled with the EP method. This means that any growth of the 
magnetic field found with the EP method cannot be due to dynamo 
action. This result is not just restricted to helical dynamos that can 
produce large-scale fields, but it also applies to nonhelical dynamos 
that produce small-scale fields. 

Major discrepancies occur even in the case of an imposed field 
and in the absence of dynamo action. It is found that the EP method 
yields excessive spectral energy, in particular at small scales. This 
discrepancy becomes even more pronounced in the dynamic case 
owing to an apparent reduction of feedback from the Lorentz force 
compared with the A method. Indeed, the saturation strength of the 
magnetic field can be about 20 times larger with the EP method 
than with the A method. 

In the ideal case, the A method may be slightly better suited to 
deal with limited numerical resolution than the EP method. How- 
ever, once the resolution becomes insufficient, there can be cases 
where, in three dimensions, spurious exponential growth occurs 
with both methods. This underlines to necessity of diffusive pro- 
cesses, but with the EP method this inevitably leads to incorrect 
solutions. 

One might expect that the A method gives good results also 
in Lagrangian schemes, because no derivative of A needs to be 
computed. An exception is the diffusion term and, of course, the 
calculation of B and J for the Lorentz force (in full MHD) and for 
diagnostic purposes. The same is true for the EP method as well. 
It should therefore be worthwhile to explore the A method also in 
Lagrangian schemes. 
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APPENDIX A: DERIVATION OF Equation © 

For completeness we give here the derivation of Equation The 
usual equation for A is 

dA 



at 



= —E — Vc 



(Al) 



where E is the electric field and (f> is the electrostatic potential. 
Using Ohm's law, J — u(E + U x B), as well as Ampere's law, 
/zn J = V x B and B — V x A we have 

BA 

— = U x V x A + ri(V 2 A - VV ■ A) - V<j>, (A2) 

where we have dropped a term (V ■ A)Vry on the RHS, because 
in our case 77 = const. Equation JA2t can be written as 

w = ~ u ^ + u ^ + ^-^- A+ ^ (A3) 

The first term on the RHS, together with the time derivative on the 
LHS, constitute the advective derivative, T)A/Dt. Next, we use 



,. dAj , dUj dUjAj 

OXi OXi OXi 

so we have 



(A4) 



DA t 1 

— = -A- ( VC7) + r,V 2 A - V(7?V ■ A U ■ A + <£).(A5) 
After a gauge transformation, A — *• A' + VA with 
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8 



A = f faV • A - U ■ A + 0) At' (A6) 
Jo 

we arrive at Equation {5}. 

APPENDIX B: DERIVATION OF EQUATION (8) 

In order to verify the R term in Equation l|9} we calculate V 2 A in 
terms of a and p, using A — | (aV/3 — /3Va), so 

+ §(a/9, w -j8a, w ). (Bl) 

Here, the last term in brackets can be written as the divergence of 
0i — ^(afljj — (3a,jj) minus h (otjPjj — P,iCx,jj) which, in turn, 
is equal to the first term in Equation l lB 1 1 , so we have 

Ai.jj = — P,jj a ,i) + ( a ,j/3,ij — P,j a ,a) + V,0i. (B2) 

The first term in brackets corresponds to the diffusion terms in 
Equation yj, the second term explains the R term in Equation 
and 4>\ gives one of several terms entering in Equation (8). 

For completeness let us here also give the derivation of the 
remaining terms. The time derivative of A is given by 

8 A 

= ±{aP,i + a/3,i - f3a,i - /3a,i) 

= 6tf3,i- 0a,i + ±S7i(af3- (3a), (B3) 

so we have 
OA 

— =dV/3-/3Va + V0 2 , (B4) 
at 

where 02 = \{°-P — P&), and dots denote partial time derivatives. 
Finally, from U ■ 'V A + A - ( VJ7) T , we have in components form 

UjAij + AjU jA = Uj(Aij - Aj,i) - X7 i( f>3, (B5) 

where 

A i,i = ~ P a ,*i) + |( Q jAi - P,ja,i), (B6) 

and 03 = U ■ A. The first term in brackets of Equation dB6l > is 
symmetric in i and j, while the second one is antisymmetric, so 
only the second one contributes to Ai t j — Aj t i, giving 

U-VA + A- (VU) T = (U ■ Va) V/3 - {U ■ V/3) Va - V0 3 . 

The first two terms explain the advection operator in Equation 10, 
while the last term contributes to — 0i + 02 + 03 in Equation {8]l. 



